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ABSTRACT 

We use high resolution Eulerian hydrodynamics simulations to study kinematic properties of the 
low ionization species in damped Lya systems at redshift z = Z. Our adaptive mesh refinement 
simulations include most key ingredients relevant for modeling neutral gas in high-column density 
absorbers: hydrodynamics, gravitational collapse, continuum radiative transfer above the hydrogen 
Lyman limit and gas chemistry, but no star formation. We model high-resolution Keck spectra with 
unsaturated low ion transitions in two Si II lines (1526 and 1808A), and compare simulated line profiles 
to the data from the SDSS DLA survey. 

We find that with increasing grid resolution the models show a trend in convergence towards the 
observed distribution of HI column densities. While in our highest resolution model we recover the 
cumulative number of DLAs per unit absorption distance, none of our models predicts DLA velocity 
widths as high as indicated by the data, suggesting that feedback from star formation might be 
important. At z = 3 a non-negligible fraction of DLAs with column densities below 10^^ cm~^ is 
caused by filamentary structures in more massive halo environments. Lower column density absorbers 
with A'^Hi < 10^^ '* cm^^ are sensitive to changes in the UV background resulting in a 10% reduction 
of the cumulative number of DLAs for twice the quasar background relative to the fiducial value, and 
nearly a 40% reduction for four times the quasar background. We find that the mass cut-off below 
which a large fraction of dwarf galaxies cannot retain gas after reionization is ^ 7 x 10^ M©, lower 
than the previous estimates. Finally, we show that models with self-shielding commonly used in the 
literature produce significantly lower DLA velocity widths than the full radiative transfer runs which 
essentially render these self-shielded models obsolete. 

Subject headings: galaxies: formation — radiative transfer — methods: numerical 



1. INTRODUCTION 

Damped Lya absorbers (DLAs) are an excellent tool 
for probing structure formation in the early Universe. 
By definition DLAs are systems with neutral hydrogen 
column densities above 2 x 10^" cm~^, and for many years 
they have been linked to forming protogalaxies at high 
redshifts. With high resolution spectrography, DLAs can 
shed light on physical processes and substructure within 
individual young galaxies. However, the exact nature of 
host absorbers has been debated for many years. 

From metal-line kinematic anal ysis of absorption 
line widths and profile asvmmetries iProc haska fc Wolfd 
p997) concluded that the observed DLA population can 
be best explained with a model in which absorption is 
caused by rapidly rotating (i^circ ~ 225 km s~^) cold thick 
{h ~ 0.3i?) disks which have already assembled at high 
redshifts, and ruled out the competing models in which 
DLAs are small (wcirc < 100 km s^"'^) protogalaxies in the 
process of hierarchical merger. However, in most models 
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of hierarchical structure formation it is extremely dif- 
ficult to have an abundant population of such massive 
disks at z ~ 3 — 4. 

On the other hand, in the high-resolution (1 kpc spa- 
tial and 5 x 10^ Mq mass resolution) SP H simulations of 
a sma ll number of high-redshift galaxies iHaehnelt et alJ 
()1998l) showed that irregular protogalactic clumps with 
I'circ ~ lOOkms"^ can reproduce the observed velocity 
width distribution and absorption profile asymmetries 
very well, without the need to invoke rapidly rotating 
massive disks. In their models large velocity widths are 
caused by a mixture of rotation, random motions, infall, 
and merging, and the resulting picture is fully consis- 
tent with most standard hierarchical structure formation 
scenarios. Despite having achieved very high numeri- 
cal resolution, this simulation features a small volume 
size (2 Mpc) and may not be representative of a typi- 
cal DLA absor ber in the cosmologica l context. Another 
shortcoming of'Ha ehnelt et all 1)1998(1 paper is that their 
models do not include radiative transfer of the ultra- 
violet background (UVB) instead assuming simple self- 
shielding above a critical density, an assumption that 
directly affects the cross-section of DLAs. In our paper 
we are going to address this issue in detail. 

Simulations with larger volume sizes inherently suffer 
from inabilit y to resolve t he l ower-mass part of the DLA 
distribution. iKatz et al.l l|199ffl studied the formation of 
DLAs in a 22.22 Mpc (comoving) volume with SPH sim- 
ulations with dark matter mass resolution 2.8 x 10^ Mq 
effectively meaning that only halos with masses above 
^ 10^^ M(7, could be resolved. 
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iGardner et all l)1997(l developed a method to account 
for absorption in halos b elow the numerical resolu tion 
of simulations. They used lPress fc Schechteil iITqtI for- 
malism to calculate the number density of unresolved 
halos and then convolved this distribution with the re- 
lation between the neutral hydrogen cross-section and 
halo circular velocity observed in their SPH simulation 
to obtain the total abundance of DLAs in a large simu- 
lation volume. We will return to the important question 
of DLA cross-sectional dependence on the dark matter 
(DM) halo velocity disper sion further in this introduc- 
tion. IGardner et al.l l|1997|) find that accounting for un- 
resolved halos with this resolution correction increases 
the DLA incident rate per unit redshift by about a fac- 
tor of two bringing it closer to the observed value but still 
underpredicting the amount of Lyman limit absorption 
b y approximately a factor of 3. 

IGardner et"airi)200H) used newer SPH simulations and 
an improved technique to compute the contribution from 
unresolved halos and found that their DLA incidence 
matches observations as long as they adopt a Vcac ~ 
60kms^^ circular velocity cutoff (corresponding to virial 
mass '--^ 2 X 10^" M0 at z = 3) below which halos cannot 
host DLA absorption systems, the value which is some- 
what higher than the more traditional threshold fcirc ~ 
40kms^^ for suppressing the formation of dwarf galax- 
ies by a UV photoionization background ijQuinn et al.l 
Il996t iThoul fc Weinbert3[T996|) . The highest dark mat- 
ter mass resolution in these simulations is 8.3 x 10^ Mq 
immediately suggesting that halos with masses below 
few X 10^" Mq are not resolved which is very similar to 
the suggested cutoff. 

On the other hand, IMaller et all l|200H) used semi- 
analytic models of galaxy formation to further investi- 
gate the hypothesis that multiple galactic clumps give 
rise to the observed kinematics and column density distri- 
butions. They studied several variants of the exponential 
disk model in which the sizes of individual galaxies are 
based on angular momentum conservation, noting that 
the resulting gaseous disks by themselves are too small 
to produce the observed kinematics and the overall num- 
ber density of DLAs. However, with a simple toy model 
they showed that a larger covering factor of the cold gas 
in these clumps can successfully reproduce the observed 
properties suggesting that lines of sight most likely go 
through tidal t ails caused by merge r s. 

Furthermore, iProchaska fc Wolfel (|2n01^ showed that 
the DLA cros s-sections ob tained from the SPH simula- 
tions (Gardn er et alJl2001j) are not consistent with the 
observed velocity width Av distribution, in fact leading 
to too many low Av systems. On one hand one would 
need a sufficient number of these low-mass systems to 
match the total observed number density of DLAs, on 
the other hand low-mass sy stems result in very small Av. 
IProchaska fc Wolf^ lj2001(l suggested two possible solu- 
tions. One possibility is that the cross-sectional depen- 
dence on the h alo mass does no t hold below the resolution 
limit of Gardner et al.l (|20nifl where processes on sub- 
kpc scales ranging from various feedback mechanisms to 
tidal stripping would disrupt this cross-sectional depen- 
dence. The second solution is that the same feedback 
mechanisms in lower-mass systems could easily double 
the typical velocity widths and are therefore central to 
understanding the overall kinematics of DLAs. 



These concerns were the focus of the recent simulations 
bv iNagamine et all l)2004D . Arguing that the relation 
between the absorption cross-section and the DM halo 
mass might not follow the same trend below 10^" M©, 
they performed a series of simulations covering a wide 
spectrum of box sizes and mass resolutions, employing a 
new conservative entropy formulation of SPH and a two- 
phase sub-resolution model for the interstellar medium 
(ISM), and varying the amount of feedback from star 
formation (SF) with a new phenomenological model for 
galactic winds. Their highest resolution run features 
216^ dark matter particles (and the same number of gas 
particles) in a 3.375 Mpc box giving 2.75 x 10^ M© mass 
resolution, whereas their largest simulation volume is 
100 Mpc on a side with 324'^ dark matter particles and 
^ 10'* times lower mass resolution. With high resolu- 
tion runs in small boxes, they see DLAs in halos with 
masses down to Mtot ^ 10^'^ Mq. Below this mass at 
z = 4.5 to z = 3 they notice a sharp drop-off in the 
DLA cross-section which they associate with photoevap- 
oration of gas in these halos by the ionizing UVB and/or 
to supernovae feedback They use the same approach 
as IGardner et alJ l|20fllf) to find the cumulative abun- 
dance of DLAs at each redshift. They fit a functional 
form to the relation between the DLA cross-section and 
the total halo mass obser ved in their simul ations, and 
then convolve it with the iSheth fc TormenI (l999) pa- 
rameterization of the dark matter halo mass function to 
obtain the total number of DLAs per unit redshift as 
a function of halo mass. They find a steeper slope for 
the DLA cros s -sectio n dependence on halo mass than 
IGardner et alJ l)2001[l resulting in fewer DLAs from low- 
mass halos, and they report a good agreement between 
the simulated and observed number of DLAs at 2; > 3. 

However, none of the simulations mentioned above ad- 
dressed both the kinematics and the gas cross-section 
properties of DLAs in the same numerical output. We 
will consider these issues in our current paper. We also 
compute self-shielding of protogalaxies from the UVB 
with a new radiative transfer algorithm. 

Before we proceed to describe our models, it is useful to 
review how the numerical simulations mentioned above 
accounted for shiel ding from the UVB bey ond the hydro- 
gen Lyman edge. iHaehnelt et alJ l)1998D have adopted 
a simple scheme to mimic the effect of self-shielding 
based on correlation between the column density and 
the physical density predicted by earlier numerical sim- 
ulations in the optically thin regime. Arguing that self- 
shielding becomes important above an Hi column den- 
sity of lO^^cm"^, which corresponds to the absorption- 
weighted density 10^'^ — 10~^ cm"'^ along the line of sight, 
they assume that all gas above a density threshold of 
10~^ cm~'^ is self-shielded and fully neutral, and lower 
d ensity gas is ionized. 

iKatz et aTl l)199(il) developed a more complex self- 
shielding correction. During actual simulation they com- 
pute the species abundances from ionization equilibrium 
in the gas that is optically thin at the Lyman limit every- 
where in the volume. Then, at the post-processing stage 
they correct projected H i maps pixel by pixel recomput- 
ing ionizational equilibrium with proper attenuation of 
the UVB in a plane-parallel slab with the same total hy- 
drogen column density and the same physical size, and 
therefore accounting for shielding in optically thick sys- 
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terns. According to lKatz et all 1)19961) . at column densi- 
ties 10^^ — 10^^ cm~^ the correction is small and varies 
between 1% and 10%, whereas at higher column densi- 
ties, A'^Hi >> 10^ ''cm"^, this correction can be as large 
as a factor of 100. Gardner ct al. (1997, 2001) all use this 
self-shielding correction analyzing models initially com- 
puted with ionizational equilibrium in an optically thin 
gas. 

On the other hand. lNaeamine et alJ 1)2004)) do not ap- 
ply any self-shielding correction at all arguing that DLAs 
with column densities above ~ lO^'' cm~^ are fully neu- 
tral even without the correction - an assumption that, in 
our view, may not hold in all situations, e.g., for dense 
shocks in tidal tails which can nevertheless be subjected 
to the photoionizing background. Instead, they use a 
more complicated two-phase ISM model, in which they 
compute the neutral hydrogen fraction with the standard 
optically thin approximation and a uniform UVB if the 
local density is below a density threshold pth which marks 
the onset of cold cloud formation. Above pth they use a 
two-phase model in which the H i fraction is a function 
of the local cooling and SF rates and a few parameters 
including feedback from supernovae. 

In this paper we present new high resolution sim- 
ulations of damped Lya systems at z = 3. To re- 
solve individual galaxies, we use a customized version 
of the adaptive mesh refinement (AMR) Eulerian hydro- 
dynamics cosmolog y code Enzo ()Brvan fc Normai3ll999t 
10 'Shea et al.ll2(l(Ml which includes simultaneous radia- 
tive transfer of the UVB beyond the hydrogen Lyman 
edge on all levels of grid refinement. Our goal is to see 
whether very high spatial resolution models with sophis- 
ticated physics can reproduce the observed kinematic and 
statistical properties of the low ionization DLA metal 
lines. 

We do not include SF feedback in the present study, 
instead focusing on the physical complications associated 
with radiative transfer in galaxy formation models with 
limited numerical resolution. We are planning to include 
feedback into our future models. 

U nlike earlier DLA simulatio ns lINaeamine 'ei~al]l200l 
iGardner et alJl2001t IHaehnelt et alJ ll998l which consid- 
ered only absorption by neutral gas in isolated (and 
mostly virialized) galaxies, we do not rule out the hy- 
pothesis that at high redshifts z > 3 a significant frac- 
tion of DLA absorption is caused by H i outside the virial 
radii of the galaxies. There is evidence of active galaxy 
assembly at those redshifts, and while DLAs are most 
likely associated with the neutral gas from which galaxies 
form, there is a certain possibility that DLA absorption 
is ca used by Hi in tidal tails in mergers fMalle r et al.l 
or even in filaments along which the gas is still 
falling into dark matter potential wells. Hachne lt et al.l 
((19985 addressed the possibility that DLAs are not nec- 
essarily caused by virialized systems, but they assumed 
a one-to-one correspondence between the local baryon 
density and the ionizational state of hydrogen which au- 
tomatically placed all neutral gas into small self-s h ielded 
prot ogalactic clumps . Sim ilarly, IGardner" eTal] (l20?nh 
and iNagamine et al.l ()2004)) identified all DM halos in 
their simulations and used a relation between the DLA 
cross-section and the total mass of associated halos to 
compute the number of DLAs per unit redshift. 

Since our models include full angular-dependent radia- 



TABLE 1 
List of models. 



Model L, h~^Mpc Base grid" A^icvcls A^grids Ax, /i^^kpc Am, Mq 

comoving at 2 = 3 comoving 
2 /i~^Mpc volume 

Al 2 128^^ 6 13,160 0.25 4.0 x lO'"^ 

4 h~-^Mpc volume 

Bl 4 128^^ 6 14,150 0.5 3.2 x 10'^ 

8 /i^^Mpc volume 

C-1 8 128^ 4 13,228 4 2.6 x 10'' 

CO 8 128=* 5 14,151 2 2.6 x lO'^ 

Cl'' 8 128^ 6 14,865 1 2.6 x lO'^ 

02'^ 8 128^ 7 15,735 0.5 2.6 x lO'^ 

"Number of base grid cells is always the same as the total number 
of dark matter particles 

''Three different UVBs were used with this model: one with the 
full stellar = 1) and full quasar (/3q = 1) components from 
Fig. (51 (CI), one with /3. = 1 and /3q = 2 (Clb), and finally one 
with /3, = 1 and /3q = 4 (Clc). 

•^In addition to the run with full UVB transfer (C2), two other 
models were computed: one with complete self-shielding above 
10~'^cm~^ (C2s), and one with the optically thin approximation 
throughout the entire simulation volume (C2t). 



tive transfer, we expect the association between the dis- 
tribution of neutral hydrogen and dark matter halos to 
emerge as we move to lower redshifts and higher levels of 
grid refinement, and consequently as more intergalactic 
HI either settles down in halos or gets destroyed by the 
UVB, but we do not put this association into our models 
in any way. In fact, the existence of high- velocity clouds 
in the halo of our own Milky Way Galaxy at present day 
(e.g. Blitz et al. 1999 and references therein) implies that 
we can expect to see extended neutral gas configurations 
at high redshifts. 

This paper is organized as follows. In section [2 we 
describe our numerical simulations, along with the two- 
stage radiative transfer algorithm, the choice of the UVB, 
and our spectrum generation routine. Sections and 0] 
present our results and conclusions. 

2. SIMULATIONS 

All simulations were run for the flat ACDM cosmol- 
ogy with an = 0.3, r^A = 0.7, f7b = 0.045, h = 0.67, 
for a primordial power spectrum (Tg = 0.9 and rig = 1. 
All models are summarized in Table 1 which lists the 
simulation volume size, the base grid size, the maximum 
number of levels of refinement, the total number of grids 
at the end of each run at z = 3, grid resolution, and mass 
resolution of each model which is simply the mass of a 
DM particle. Most of the models include full radiative 
transfer, except for runs C2s and C2t which were com- 
puted with self-shielding above density 10~^ cm~^ and in 
the optically thin regime, respectively. 

We identify the following three parameters in our cal- 
culations: grid resolution, box size/mass resolution, and 
the UVB amplitude. We consider the box sizes of 2h~^ , 
Ah~^ and 8/i~^ Mpc, with base grid resolution 128'^ and 
up to seven levels of refinement by a factor of two. In 
Enzo the total number of DM particles is always the 
same as the number of base grid cells N^^^ resulting in 
mass resolution 
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AM « 3.2 X 10^ f f^^^ 'mq. (1) 

Therefore, for a fixed base grid, mass resolution is deter- 
mined by the current box size. 

2.1. Two-stage radiative transfer 

We use a two-step algorithm for radiative transfer. 
First, we evolve all simulation volumes solving the radia- 
tive transfer equation in a low angular resolution mode 
simultaneously with the equations of hydrodynamics. To 
correct for low resolution, we apply a high angular reso- 
lution transfer filter iteratively to our solutions at z = 3 
to improve on the equilibrium position of HI, Hel and 
HeH ionization fronts. 

We developed a UVB radiative transfer module for the 
parallel AMR cosmology N-body Eulerian hydrodynam- 
ics code Enzo. The fluid flow equations are solved with 
the PPM scheme ( ColeUa & Woodward 1984) on a co- 
moving grid (see for iO'Shea et al.l l|200 4'l for details), and 
chemical abundances for various ionization and molec- 
ular states of hydrogen and helium are solved with a 
chemi cal reaction network with 9 species a nd 28 reac- 
tions llAnninos et alJll997t lAbel et alJll997D . Each level 
of grid refinement has twice the spatial resolution of the 
previous level. We solve the radiative transfer equation 
with a photon conserving scheme self-consistently at each 
level of resolution carrying fluxes explicitly from parent 
grids to all of its subgrids. On each grid, radiative trans- 
fer is computed with a timestep which is usually signif- 
icantly smaller than the hydro timestep and is adjusted 
adaptively to obtain the best balance between accuracy 
and the speed of calculations. Due to complexity of the 
setup we have not included point source radiation in our 
models yet, and also we limit transport of background ra- 
diation to simple sweeps along the xyz-coordinate axis. 
We then post-process Enzo output with much higher an- 
gular resolution transfer. Below we describe these two 
steps in detail. 

2.2. Stage one: time stepping and parallelization of the 

base grid 

By definition, radiation propagates at the speed 
of light. In an explicit advection numerical scheme 
the Courant condition would necessarily require pro- 
hibitively small timesteps to guarantee stability. How- 
ever, our photon conservation technique is inherently sta- 
ble, and from this standpoint there is no need to take very 
small timesteps, but accuracy of the solution is an en- 
tirely different issue. In many cases where the supply of 
ionizing photons greatly exceeds the recombination rate 
I-fronts can propagate at or close to the speed of light, 
and even in this regime our photon conservation will give 
an accurate solution as long as the radiation-chemistry 
timestep is small enough. How do we ensure that we 
provide timesteps small enough to balance accuracy and 
speed of calculations? 

Let us first describe our numerical setup. Since in gen- 
eral radiation driven fronts can propagate much faster 
than the fluid flow, radiative transfer and chemistry nor- 
mally have to be iterated multiple times per hydro time 
step. On each subgrid, independently of the resolution 



level, first we perform a hydrodynamical update at a 
fixed timestep Atf determined from the hydrodynami- 
cal Courant condition on that subgrid. We then proceed 
to compute photon conservation in all directions on a 
much smaller timestep Aij. — Att/(2" — 1) where n is 
some positive integer number. We then update temper- 
ature and solve the chemical rate equations on the same 
small Atf. Next we compare distributions of some state 
variable, e.g., a fraction xmi of ionized hydrogen on our 
subgrid. If the maximum change in xmi during Afi- on 
the subgrid does not exceed 10% then we increase Afr 
by a factor of two, otherwise keep the same small Atj-, 
and do another radiation-chemistry update, and repeat 
the whole procedure of these updates until we reach the 
end of the hydro step Atf . Ideally, if there are no I-fronts 
on the subgrid we can do the entire radiation-chemistry 
calculation with n subcycles for each hydro step. On the 
other hand, in the worst case scenario with fast I-fronts 
we will need (2" — 1) subcycles. The exact number of 
subcycles between n and (2" — 1) is determined automat- 
ically by the code based on maximum relative changes in 
xmi in our subvolume. 

To determine a suitable value of n, we ran a number of 
tests with Air = Aif/(2" — 1) fixed throughout the entire 
run from high to moderate redshifts. We observe fast 
convergence as the number of subcycles reaches few tens, 
but the solution is accurate enough even for (2" — 1) ^ 10. 
In all our calculations we use n = 4 varying the number 
of subcycles between n and (2" — 1) as described above. 

Ionizing photons in our model all originate at the edge 
of the base grid, and propagate inward. Parallelization 
in Enzo is done through volume decomposition in which 
the lowest resolution base grid is subdivided equally be- 
tween all processors. Since we treat radiation explic- 
itly, solutions to the radiation field in all directions on 
all processors have to be connected to each other at ev- 
ery radiation-chemistry subcycle, and at the same time 
we need to minimize the amount of idle time each pro- 
cessor spends waiting for flux updates from neighbor- 
ing volumes. To achieve this goal of optimal paralleliza- 
tion and best load balancing while simultaneously trans- 
porting photons in all directions, we constructed the fol- 
lowing algorithm. At the beginning of every radiation- 
chemistry subcycle each processor scans all angular direc- 
tions in some preset order checking if input fluxes (from 
the boundary or from another processor) are available in 
that direction and if radiative update has not been done 
yet. When it finds such an angle, it performs radiative 
transfer in that direction, while all other processors are 
doing their own updates. If no input flux in any direction 
is available then the processor sits idle for the duration 
of this directional subcycle. When the subcycle is done, 
fluxes are passed to neighboring volumes, and a new sub- 
cycle begins. To further illustrate this idea, we drew a 
sequence of directional updates for the 3x2x1 volume 
decomposition in Fig. ^ In this particular example we 
have a 100% parallelization efficiency, which is clearly 
not always achievable for larger decompositions. 

To summarize, in the hierarchy of events on a base grid 
every single hydrodynamical step is followed by a series 
of radiation-chemistry updates each of which in turn con- 
tains directional subcycling with flux exchanges among 
individual processors. When a spatial region is reflned, 
radiation fluxes for both stellar and quasar components 
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Fig. 1. — A sequence of directional updates for the 3x2x1 
volume decomposition. 




in all six directions are passed from parent grids to sub- 
grids, along with all hydrodynamical variables. 

2.3. Stage two: high angular resolution transfer 

Our radiation hydrodynamics (RHD) calculation is 
parallelized via domain decomposition, in which we pass 
radiative intensities from one processor to another on 
the base grid, and in turn on each processor we advect 
these intensities along the AMR grid hierarchy. This ap- 
proach results in a fairly complex algorithm, forcing us 
to limit RHD transport only to the six directions along 
the major axes. The total number of photons entering 
the volume in plane-parallel fluxes is kept the same as if 
we had an isotropic incident distribution, therefore those 
photons which burn their way into a denser environment 
at small solid angles tend to over-heat and over-ionize 
the medium. At the same time this approximation nat- 
urally leads to formation of shadows behind self-shielded 
regions which in reality could be exposed to ionizing ra- 
diation. To compensate for this shortcoming, wc post- 
process our 3D ionization and temperature distributions 
with a much higher angular resolution radiative trans- 
fer code based on a n ew fully threaded tran s port e ngine 
algorithm (FTTE) bv'Razo umov fc Cardall l|2005j) . 

At each output redshift we convert the block- 
structured AMR output of our coupled RHD runs to a 
fully threaded format and use it as a first guess to com- 
pute radiative transfer with FTTE along 192 directions 
chosen to subdivide the entire sphere into equal solid an- 
gle elements. We then use this updated radiation field 
to compute the ionizational structure and temperature, 
and then iterate in radiative transfer and chemistry un- 
til we find the equilibrium positions of HI, Hel and Hell 
ionization fronts at a given redshift. We find that ~ 20 
iterations are needed for convergence. 

2.4. Adopted UVB and frequency dependence 

Constraints on the UV B at high redshifts come from 
the Lya forest studies. iScott et all 1)20001) used the 
proximity effect in quasar absorption spectra to de- 
rive the UVB amplitude assuming that it has the same 
spectrum as individual quasars. Using Ly-alpha emis- 
sion line redshifts, they get the value Ji, = 1.4to'5 x 
10-21 erg s-^cm-^Hz-^sr-i for 1.7 < z < 3.8. On 
the other hand, with OIII and Mgll redshifts they get a 
lower value J„ = 7.0144 ^ lO^^^ ergs^^ cm^^ Hz"^ sr"^ 
for 1.7 < z < 3.8, corresponding to HI photoionization 
rate 1.9t}:^ x 10" ^^ 8"^ 

More recently, iTvtler et all 1)20041) and iJena et all 



Fig. 2. — Amplitude of the assumed fiducial UVB at 13.6 eV 
in standard units of 10"'^^ ergs~^ cm~-^ sr~^ Hz~i. Shown are the 
quasar (dashed line) component /2i,q and the stellar (dotted line) 
component /21,*, and their sum (solid line) at the hydrogen Lyman 
edge for ff * = = 1. For compa rison we also plotted the total 
UVB from lAbel fc HaehnelH IT999h (dash-dotted line). 

()2005)) build a concordance model of the Lya forest at 
z = 1.95 using a series of hydrodynamical simulations 
with grid sizes up to 1024"^ and box sizes up to 76.8 Mpc 
to reproduce the observed flux decrement from the low- 
density intergalactic medium (IGM) alone. Their UVB 
corresponds to an ionization rate per Hi atom of r9i2 = 
(1.44 ± 0.11) X lO-^^s"^, which is slightly high er than 
the earlier estimates, e.g., the prediction bv lMadau et alJ 
(Q.,999) with 61% from QSOs and 39% from stars. The 
re dshift evolution of th i s conc ordance model can be found 
in Paschos & Norman' (' 2005D . A more recent study by 
Kirkman et al. (2005) extends this higher estimate of 
r9i2 ~ 1.4x lO^i^ s^i to the redshift range 2.2 < z < 3.2 
translating into ~ 5 x 10-^^ ergs"^ cm-^ Hz^^ sr^^. 

The background in the vicinity of DLAs in a cluster 
that does not host a quasar is most likely going to be 
smaller than these average values, due to the geometrical 
dilution factor, and local attenuation in the cluster. In 
our other models not listed in Table 1 we experimented 
with a variety of UVB spectra and amplitudes, and as a 
reference we decided to use a smaller background based 
on Abel & Haehnelt (1999). We included both stellar 
and quasar ionizing photons adopting a two-component 
power-law UVB 

L = (j^^ + /5q/21,q (j^^ ' (2) 

where a* = 5, aq = 1.8, and /3* and /3q are the model pa- 
rameters. We extend the stellar and quasar components 
of lAbel fc Haehnelti l)1999D to higher redshifts with the 
expressions 

/2i,q = Q + a(^-Q), (3) 
hi,* = (5 + a{Q - S)) (1 - tanh(0.5z - 7)) /2, 
where 

10exp[-(z/2.5)3] 
^ l + (7/(l + z))4 ' 
1-b exph(z/4)3] 
2 l + (7/(l + z))4 ^ 
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b (l + z)3-3 5 
106.4 
a = 0.3 exp[— (z 



- exp 



-(z-0.5)^ 



l + (z + 2.09)2075/16 



4.5)V4], 
6=l + tanh(1.5z-6.3). 



We plotted /21,*, ^21, q and the background from 
Ebel & Haehnclt' (1999)' in Fig. [3 We allow for a sub- 
stantial stellar component at z > 8, which is consis- 
tent with the current theoretical v iews of the global his- 
tory of SF in a ACDM Universe ijSpringel &: Hernauisti 
|2003^ in light of the discovery of a large optical depth 
to electron scattering by WMAP (Kogut et al. 2003J 
suggesti ng presence of bright ionizing sources at z ~ 
15 - 17 llWvithe fc LoeHmM IHaiman fc Holder! I2?M 
iHui fc Haimanll2003t iCiardi et al.ll2003j) . As the main 
focus of this paper is to study the resolution effects of 
radiative transfer in galaxy formation models, we do not 
explore other possible star formation histories, but un- 
doubtedly the cumulative history of star formation at 
2 > 6 might have an important effect on the observable 
properties of galaxies even at lower redshifts. 

We use three frequency bands, 13.6 — 24.6 eV, 24.6 — 
54.4 eV, and above 54.4 eV, to compute Hi , Hei and 
Hen ionization. Inside each frequency band [vgyVg+i] 
the transport variable is the intensity at frequency Vg 
weighted by a corresponding power-index-dependent fac- 
tor 



1 - (Vg+l/Vg)^ 

a. - 1 



(5) 



where g — 1,2,3. All radiation-related rates (photoion- 
ization, photoheating, and absorption) also depend on 
Qfg. Finally, inside each frequency band we use a sin- 
gle transport variable describing both stellar and quasar 
components with an effective index Ug constructed to en- 
sure that the total energy inside the band is equal to the 
sum of the two individual components. 

2.5. Spectrum generation 

To analyze the statistical properties of DLAs, we drew 
100,000 random (through a random point in a random 
direction) lines of sight through the entire grid hierarchy 
of each model using the highest resolution cells available 
at each point along the line, and studied all absorption 
systems with HI column densities iVni > 2 x IQ^'^ cxnT^ . 
For analysis we use two "low ion" lines of Si II at 1526A 
and 1808A. Since the latter line has a significantly 
lower oscillator strength, we use it only in absorbers 
with log(A'Hi/ cm~2) > 20.6, and for the lower column 
densities use the former transition. We assume a uni- 
form abundanc e throughout the volume with metallicity 
[Si/H] = -1.3 llWolfe et al.ll2005D . Noise is added to our 
artifical spectra such that the resultant signal-to-noise 
ratio of each Ikms^^ pixel is 20:1. We define the Si II 
line width as the width of the central part of the profile 
responsible for 90% of the integrated optical depth as 
described in lProchaska fc Wolfd l|1997j) . In the unlikely 
case that a line of sight crosses multiple absorbers, we 
consider two components to be caused by a single DLA 
if the separation between the centers of the two corre- 
sponding line profiles is less than 400kms~^, otherwise 
we just analyze the strongest component. 



3. RESULTS 

To demonstrate the nature of DLAs produced in our 
models, in Fig. Owe plotted the projected HI column den- 
sity in a volume 8 Mpc on a side and 8 Mpc thick, for run 
CI at z = 3, with zoom-ins on a cluster of galaxies and a 
disk galaxy. All objects with HI column densities above 
2 X 10^" cni^2 give rise to DLAs. To investigate the va- 
lidity of our results and compare them to observations, 
we use two types of distributions: the HI column den- 
sity frequency distribution f{N,X) and the hne density 
^dla(-^) of DLAs with the Si II velocity width higher 
than vsiu vs. vsai- The frequency distribution f{N,X) 
is defined such that f{N, X)dNdX is simply the number 
of DLAs in the intervals [N, N + dN] and [X, X + dX], 
where dX is the "absorption distance" interval 



dX = ^" 



H{z) 



(1 + zydz. 



(6) 



In principle, any combination of the following four fac- 
tors can affect our results: finite grid resolution, finite 
mass resolution, the amplitude and spectrum of the as- 
sumed UV background, and radiative and mechanical 
feedback from SF. As the goal of this work is to demon- 
strate a method to compute the ionization structure of 
the outskirts of high redshift galaxies with self-consistent 
radiative transfer of the UVB, we do not include the un- 
certainties of feedback from SF here. In this paper we 
concentrate on resolution issues in the new models and 
the dependency on the UVB, at a fixed redshift (z = 3). 
In our next paper we will study the redshift evolution 
and stellar feedback. 

3.1. Numerical resolution 

In Fig. 01 we plotted f{N,X) for all of our runs at 
z = 3. For comparison, we also plotted the data from 
the SDSS DLA survey from iProchaska et alJ (|2005H for 
redshift intervals z = 2.5 — 3.0 and z = 3.0 — 3.5. The 
top panel shows our 8h~-^ Mpc models, for which the grid 
resolution ranges from 4/i~^ kpc (comoving) for run C-1 
to 0.5/i~^ kpc for run C2. At this redshift these numbers 
translate into 1.5 kpc and 190 pc physical resolution, 
respectively. Although the results have not converged 
yet at our highest resolution, the trend is clear: as we 
increase grid resolution, the baryons in halos can col- 
lapse further yielding smaller DLA cross-sections. For 
any curve, the disagreement between models and obser- 
vations is the largest for high-column density systems 
which have the sharpest baryon concentrations and are 
particularly sensitive to grid resolution. 

The challenging aspect of the DLA modeling is clearly 
illustrated in the lower panel in Fig. 01 in which we apply 
the same 128'^ base grid with 6 levels of refinement to a 
smaller (4/i~^ Mpc) volume, resulting in two times better 
grid resolution and eight times higher mass resolution 
(run Bl, dashed line). Unlike model CI (solid line), this 
run has a population of self-shielded halos in the 10* — 
10^ Mq mass range. While naively one would expect to 
see higher f{N,X), especially at low column densities, 
this eiTect is almost precisely offset by the reduced size 
of individual halos through higher grid resolution. This 
is further seen in run Al (dotted line), which extends the 
population of self-shielded halos down to ~ 7 x lO*" Mq 
but features almost identical f{N,X). 
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Fig. 3. — HI column density in a volume 8h ^ Mpc on a side and 8h ^ Mpc thick (comoving), from run CI at z = 3. The zoom-in boxes 
are 600?t~^ kpc (top) and 80h~^ kpc (bottom) on a side in comoving units, translating into an approximately lOkpc (physical units) DLA 
around the disk galaxy in the center of the bottom panel. The DLA column densities correspond to green and up in the color map on the 
right. The resolved dynamic range in this model is 8192. 



The same effects can be observed in Fig.^lwhich shows 
the Une density d-DhAiX) of DLAs per unit absorption 
distance with the Si II velocity width higher than vsai- 
This figure is similar to the cumulative a bundance plots 
vs. h alo ci rcular velocity ( o r mas s) in IGardner et al.l 
1)20011) and iNaeamine et aLl l|2004|) . but we chose the 
neutral gas line width as a more generic measure of the 
strength of the absorber, since our simulated absorbers 
can also be caused by tidal tails and filaments. The hor- 
izontal lines in Fig. [S] s how th e me an observed DLA line 
density from IProchask a et al.l 1)200 5.') . 

In the top panel we see a monotonic decrease in the 
line density of DLAs as we increase the grid resolution. 
The higher velocity tail is not affected until we compute 
our highest resolution model C2 which starts to resolve 
individual minihalos in the highest density regions. With 
a certain probability that is related to the physical size of 
their neutral cores, these small galaxies (with a few tidal 
streams mixed in) can be seen as multiple absorption 
components in a single line of sight to a remote quasar. 

This effect can be further illustrated in Fig.Elin which 
we plotted Si II 1526^ or 1808^ line profiles for a typical 
DLA (top profile) and for nine DLAs with largest velocity 
widths (profiles 2-10, from top to bottom) in model C2. 
With the line detection criteria defined in Sec. 12.51 we 
found the following velocity widths for these spectra: 31, 
468, 191, 149, 281, 168, 170, 306, 193, and 205 kms'^ 



(from top to bottom panels in Fig. |SJ). Spectra 2 and 
8 (counting from the top) are particularly clear exam- 
ples of several components falling onto the same line of 
sight. These multiple component DLAs give rise to dis- 
tinct tails at higher velocities which stand out in our high 
(190 pc physical) resolution models C2 (upper panel) and 
Bl (lower panel). While the highest resolution run Al 
(lower panel) shows a similar tail, it is shifted towards 
lower (65 — lOOkms"^) velocity widths, in part due to 
the smaller velocity dispersion in the cluster, and in part 
due to the lower cross-sections of individual absorbers. 

To summarize, none of our models at this stage are 
able to reproduce the observed velocity width distribu- 
tion for systems with vsni ^ 30kms~^. Although the 
resolution effects are complex, and none of our models 
have fully converged to the observed column density dis- 
tribution yet, it seems plausible that some of the effects 
described in this section can drive the line profile dis- 
tribution to even lower velocities. For example, as we 
increase grid resolution, fewer and fewer systems will be 
crossed by the same line of sight, moving many DLAs 
in Fig.Elfrom 100 -400 km s^^ into the 30- 100 km s'^ 
range. Wc therefore conclude that feedback from star 
formation seems to be the most likely mechanism to get 
DLA velocity widths as high as indicated by observa- 
tions. In our future simulations, in addition to exploring 
feedback, we will also strive to increase mass resolution 
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Fig. 4. — Top panel: the upper four solid lines show the HI 
column density frequency distribution f{N, X) for runs C-1 , CO, 
CI, and C2 (from thin to thickest line), all at 2 = 3. The lower 
two lines are the F-func tion fits to the SDSS DLA survey data from 
IProchaska et al.l ^2005^ for z = 2.5-3.0 (dashed) and z = 3.0-3.5 
(dotted). Lower panel: same as in the top panel, except that the 
upper three lines give f{N,X) for runs CI (solid), Bl (dashed), 
and Al (dotted). 



Fig. 5. — Top panel: line density lohAiX) of DLAs with the Si 
II velocity width higher than vgm vs. fsiii, at z = 3. The lines 
correspond to models C-1 , CO, CI, and C2 (from thin to thickest 
line). The dotted line shows the observed kinematic distr ibution 
from d ata at all redshifts, compiled from Fig. 10 in IWolfe et alj 
120051) and normalized to match the cumulative observed value 
0.067 ± 0.006 at z = 3 (horizontal line). Lower panel: same as in 
the top panel, except that the three thick lines give ^dla(^) for 
runs CI (solid), Bl (dashed), and Al (dotted). 



in progressively larger simulation volumes, which should 
produce many more self-shielded halos with masses of 
few 10® Mq in galaxy clusters with larger velocity dis- 
persions. 

3.2. Halo and intergalactic DLAs 

We find that the nature of DLA absorption is a func- 
tion of the halo environment. The cumulative halo mass 
function in our 8h~^^, Ah^^ and 2/i~^Mpc volumes at 
z = 3 is plotted in Fig. [3 Most of these halos give rise 
to DLAs; however, not all DLAs can be associated with 
individual halos. 

Gardner ct a^ ljl997j) assumed a relation between the 
neutral hydrogen cross-section and the host halo mass to 
predict the total abundance of DLAs in a large simula- 
tion volume. On the other hand. iNaeamine et al.l l)2004f) 
computed this relation for individual halos in a series 
of SPH models with the assumption that all DLAs are 
fully neutral and therefore a self-shielding correction is 
not necessary. 

Since we include the effects of the UVB on gas at the 
surface of DLAs, it is interesting to make an independent 
estimate of neutral hydrogen cross-sections for individ- 
ual absorbers in our models. We draw a large number 
(A'los — 10®) of random lines of sight parallel to all three 
major axes, and find all absorption systems with the neu- 



tral hydrogen column density above 2 x 10^°cm~^. We 
then look for all halos inside the same base grid cell as 
the geometrical center of a DLA, find the one which is 
closest to this DLA, and associate it with this absorption 
system. Some extended DLAs caused by intergalactic HI 
will have no hos t halos. Halos were ide ntified with the 
HOP algorithm ijEisenstein &: Hul]ll998j) using the rou- 
tine enzohop which is part of the Enzo package. Follow- 
ing this procedure, we count the number of damped sys- 
tems iVDLA associated with each halo. Then the effective 
radius of the absorber can be estimated approximately 
as 

r.. = a^/-'^L,^Jp^Y\ (7) 

V ^VLOS / 

where Lbox is the physical box size. In Fig.|Slwe plotted a 
relation between the halo mass Mhaio and its absorption 
radius rcff in runs Al (diamonds), Bl (triangles), and C2 
(crosses) . A least-squares fit to these data with equation 

log(reff/ kpc) = a log(A/haio/ Mq)+ 13 (8) 

gives a ~ 0.38 and ~ —3.07 (thick solid line). At 
z — 3 our absorption cross-sections ar e approximately 
a factor of two lower than the ones in iNagamine et alJ 
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Fig. 6.— Si II 1526A (or 1808A) line profiles for ten selected 
DLAs in run C2 at z = 3. 
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Fig. 7. — Cumulative halo mass function at z = 3 for dark 
matter only, for runs CI (thick solid line), Bl (thick dashed line), 
and Al (thick dott ed line). The thin solid line is the analytical 
mass function from lMo fc Whitd 120021) . 



l)2004(l . To avoid any confusion, we want to stress that 
this finding does not contradict our DLA counts in Fig.^ 
and 13 since the halos in Fig.|Sldo not represent all DLAs 
in the volume as this plot does not include intergalactic 
HI clouds. In addition, our data in Fig. |H1 should be 
viewed as lower limits to DLA cross-sections rather than 
their actual values since in a massive cluster environment 
very often we cannot identify a single halo responsible for 
a particular DLA. For example, in Fig.|31we can see larger 
HI clouds engulfing multiple halos. For this reason, in 
our estimate of r^s we only searched for halos inside the 
same base grid cell as the geometrical center of a DLA 
which undoubtedly underestimates the true absorption 
cross-section. 

In addition, interactions between galaxies in more mas- 
sive environments create tidal streams which can pro- 
duce DLA lines. The cluster in our 8h~^ Mpc vol- 
ume at z = 3 has extended HI tails formed by galaxy- 
galaxy interactions which - at certain line of sight ori- 
entations - produce column densities in the range 2 x 
10^° — 5 X lO^-^cm"^. The typical physical density in 
the self-shielded filaments in these streams is of order 
~ 10~^cm~'^, the hydrogen neutral fraction often (but 
not always) exceeds 0.9, and physical widths are of order 
5 — 10 kpc. Gas which is not self-shielded from the UVB 
will never reach these column densities, unless perhaps 
it is highly compressed into shocks by feedback which we 
do not include here. 

To estimate the fraction of the line density ^dlaI-'^) 
due to gas outside of halos, we removed all neutral hy- 
drogen from inside the virial radii of all halos in the high- 
est resolution 8h~^ run C2, and reran the analysis. We 
defined the virial radius as the radius of the sphere that 
encloses 180 times the mean mass density of the Universe 
at that redshift. We found that at our highest numerical 
resolution approximately 29% of all DLAs are due to gas 
outside of galaxies at 2: = 3. While in the original run 
C2 the highest DLA column density is lO^^-^cm"^, it 
changes to 10^^'^ cm~^ if we consider only intergalactic 
gas. 

On the other hand, in less massive halo environments, 
especially in the smaller 4ft,~^ and 2/i^^Mpc simulation 
volumes, the main reason behind our fairly large DLA 
counts despite the small effective cross-sections is a much 
larger contribution to f{N,X) from lower mass halos in 
the range ~ 7 x lO'' — 10^ M© which we demonstrate 
below. 

3.3. Photoionization of low-mass galaxies 

It is well known that low-mass galaxies cannot retain 
gas after reionization. The exact value of the cut-off is 
not very well establishe d and without doubt de pends on 
the local environment. iNaeamine et alJ l|2004fl notice a 
sharp drop-off in the DLA cross-section below the mass 
~ IQ^Mq which they associate with photoevaporation 
of gas in these halos by the ionizing UVB and/or with 
supernovae feedback. 

Among all of our models, the least massive halo which 
is associated with a DLA has a dark matter mass 3.6 x 
10^ M© (Fig. It was found in run Al in which the 
halo mass function extends down to 1.4 x 10^ M© (Fig.EJ 
which is ~ 35 times the mass resolution of this run. The 
two other least massive halos associated with DLAs were 
also found in this highest resolution run. However, it is 
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Fig. 8. — Effective radius vs. halo mass, for all DLAs with 
identified host halos in runs Al (diamonds), Bl (triangles) , and 
C2 (crosses) at z = 3, for the fiducial background in Fig. This 
figure does not account for DLAs caused by intergalactic HI clouds 
for which no host halo in the same base grid cell was found. The 
dotted horizontal line shows the size of the base grid cell in the 
8h~^ Mpc volume. The solid lines are our least-squares fits with 
eq. |S]for all halos from runs Al, Bl and C2 (thick line), and all 
halos from runs Al, Bl and CI (thin line). For comparison, we 
also plotted the fits to DLA radii from Nagaminc et al. (200'f), for 
models 03 (dashed line, no wind, 10 Mpc volume), Q3 (dashed- 
dotted line, strong wind, 10 Mpc volume, low resolution) and Q5 
(dotted line, strong wind, 10 Mpc volume, high resolution). 



only above the mass ~ 7 x 10^ Mq that we start seeing 
a significant population of halos with HI in absorption. 

The fact that we see neutral gas i n halos with m asses 
slightly below the cut-off of Nagami ne et al.l l|2004f) can 
be attributed to both the spatial variations in the UVB 
inside the galaxy cluster and the lack of supern ova feed- 
back in our models. As iNaeamine et al] l)2004j) pointed 
out, a sharp cut-off should be expected due to the strong 
dependence of baryon cooling and photoheating on the 
virial temperatures of halos around ~ lO^K. In fact, 
we see such a cutoff at '--^ 7 x 10'' Mq in Fig. |S1 As all 
photoionizing radiation in our models comes from the 
box boundary, few low-mass halos might be shadowed 
by other more massive systems and be able to retain 
their neutral gas, as we see in two or three galaxies with 
masses below 7 x 10^ M0. Furthermore, a number of 
galaxies in the mass rang e 7 X 10'^ Mq - 10^ Mq are not 
exposed to the full UVB as part of their "sky" is blocked 
by nearby absorption systems (see, e.g., the zoom-in pan- 
els in Fig. |2Jl, and the mass cut-off above which galaxies 
can retain neutral gas is shifted towards slightly lower 
masses. Of course, this effect depends on the environ- 
ment and will be more pronounced at higher redshifts 
where the spatial variations in the UVB are larger. Star 
formation and supernova feedback will partially reverse 
this effect removing neutral material from lower mass 
galaxies assuming that these galaxies have accumulated 
enough cold gas to host star formation in the first place. 
However, this effect is very complex and, in addition to 
the mechanical feedback from winds and supernovae, will 
have to include transfer of locally generated UV photons. 

3.4. Dependence on the UVB 

Our distributions are somewhat sensitive to a change 
in the UVB. In model CI we varied the quasar compo- 
nent from the fiducial value in Fig. |21to twice (Clb) and 
four times (Clc) that amplitude and plotted the results 
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Fig. 9. — The three thick lines in each panel show results for a 
run with the full quasar background from Fig. |^ (CI, solid line), 
and runs with two times (Clb, dashed line) and four times (Clc, 
dotted line) the fiducial quasar background, all at z = 3. Top 
panel: HI column density frequency distribution f{N,X). The 
two thin lines are the SDSS DLA survey data as in Fig. |4] Lower 
panel: line density ^dla(^) of DLAs with the Si II velocity width 
higher than vsm vs. vsm. The thin dotted and thin horizontal 
lines show the observed kinematic distribution as in Fig. |S] 



in Fig. |5| As we increase the UVB amplitude and si- 
multaneously steepen its spectrum, more HI is ionized 
everywhere on the outskirts of galaxies, but in terms of 
DLA statistics - not surprisingly - primarily low col- 
umn density systems are affected. While the under- 
resolved model CI features the cumulative number of 
DLAs ~ 25% above the observed value at z = 3, the 
higher UVB model Clc has the number of DLAs ~ 25% 
below the observed value. Noticeable differences between 
the two models extend to the column density 10^^'^ cm~^, 
but the higher column density systems are not affected. 



3.5. Effects of radiative transfer 

In addition to the full radiative transfer run C2 in 
the 8ft.^^Mpc volume, we also computed a model C2s 
with complete self-shielding above the physical density 
10^^ cm^^, assuming that the local UVB is just the uni- 
form cosmic background from Fig. [21 in every cell be- 
low this threshold, and zero above it. This assump- 
ti on has been previously used in DLA modehng, e.g., 
in iHaehnelt et al.l l)1998|) . We also computed a model 
C2t with the optically thin approximation imposing the 
same uniform cosmic UVB in every cell in the volume. 
In Fig. ^1 we plotted distributions for these three mod- 
els. All three models were iterated until H/He ionization 
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equilibrium. 

It is evident that self-shielding of neutral gas is 
the dominant mechanism determining all properties of 
DLAs. The model with the uniform background in the 
optically thin regime produces very few DLAs as ex- 
pected. Surprisingly, the full transfer model and the 
one with self-shielding give roughly the same cumulative 
number of DLAs. However, since the transfer model in- 
cludes high energy photons, the ionizational structure of 
HI regions is more complex than in the fully shielded 
case. These differences result in vastly different line 
widths distributions: in the self-shielded model most 
lines widths are clustered around 15kms~^, whereas in 
the model with transfer they are concentrated at higher 
velocities in the range 15 — 30kms~^. It is reassuring 
that the model with full transfer in the lower panel in 
Fig-Elis about half-way from the self-shielded model to 
observational data. 

It is important to point out that our self -shielded model 
produ ces results very different from iHaehnelt et al.l 
l)1998j) . namely substantially smaller velocity widths on 
average. A number of factors could contribute to this dif- 
ference, such as the use of Eulerian hydrodynamics with 
AMR instead of SPH, but perhaps the most important 
factors are our much larger sample of galaxies (we used 
710 distinct halos to generate 406 DLAs in runC2s vs. 
their sample of 40 protogalactic clumps used repeatedly 
to produce 640 DLAs), and our much higher grid res- 
olution (190 pc physical vs. their 1 kpc). But as we 
point out, our full radiative transfer runs produce fur- 
ther improvement in kinematics modeling rendering any 
self-shielded models obsolete. 

In Fig. ^2 we plotted the mean specific intensity at 
13.6 eV, 24.6 eV and 54.4 eV vs. the local gas number 
density at z = 3. One can easily see hardening of the 
spectrum in self-shielded regions above the physical den- 
sity lO^^cm"*^, as more energetic photons travel further 
in the neutral medium. 
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Fig. 10. — The three thick lines in each panel show results 
for a run with full radiative transfer (C2, solid line), a run with 
self-shielding above the physical density 10~^ cm~^ (C2s, dashed 
line), and a run with the optically thin approximation throughout 
the entire simulation volume (C2t, dotted line), all at 2: = 3. Top 
panel: HI column density frequency distribution f(N, X). The two 
thin lines are the SDSS DLA survey data as in Fig.|4] Lower panel: 
line density iwLAi^) of DLAs with the Si II velocity width higher 
than vsiii vs. usiii- The thin dotted and thin horizontal lines show 
the observed kinematic distribution as in Fig. ISl 



4. SUMMARY 

The standard hierarchical model is very successful in 
explaining the emergence of structures in the early Uni- 
verse. DLAs are viewed as neutral gas clouds confined 
to individual galaxies, whether in virialized systems or 
systems still in the process of merging. This paradigm is 
clearly seen in most previous DLA simulation literature 
where analysis is based on a fit to the relation between 
DLA cross-section and h alo mass. The results of high 
resolution simulations bv IHaehnelt et all l)1998|) suggest 
that small protogalactic clumps in the process of merger 
could explain the observed line profiles, and the statis- 
tical analysis of abso rption line kinematic s within semi- 
analytical models bv iMaller et al.l l)2001tl suggests that 
tidal tails can be responsible for absorption. 

To test both these hypotheses and a broader idea that 
DLAs can be caused by absorption in neutral intergalac- 
tic clouds either stripped of galaxies during interactions 
or even still accreting onto galaxies from the IGM, we 
propose to include radiative transfer of the UVB into 
galaxy formation models, and argue that this new piece 
of physics is essential in order to reproduce DLA ob- 
servables such as the column density and velocity width 
distributions. An accurate treatment of background pho- 
tons allows us to address the ionization structure of the 
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bg total gas density, cm"^ 

Fig. 11. — Mean specific intensity (ergs~^ cm~^ Hz~^ sr~^) vs. 
gas number density in model CI at ^ = 3 at photon energies 13.6 eV 
(solid line), 24.6 eV (dashed Une) and 54.4 eV (dotted line). 



outskirts of high redshift galaxies in a much more precise 
way. 

We simulate galaxy formation in a series of cosmolog- 
ical volumes ranging in size from 2 Mpc to 8 Mpc, solv- 
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ing simultaneously for the first time coupled equations 
of hydrodynamics, self-gravity, multi-species chemistry 
and radiative transfer to account for shielding against 
the UVB. Radiative transfer is computed on all levels 
of grid refinement with a two-stage calculation, first just 
along the three major axes simultaneously with hydro- 
dynamics, and then in 192 directions in the postprocess- 
ing mode. We use AMR to resolve individual galaxies 
with the same r efinement criteria every where in the vol- 
ume. Similar to'Naga mine et aP l)2004f) . we use a series 
of computational volumes of different sizes while varying 
numerical resolution to try to achieve optimal balance 
between resolution and fair volume sampling. However, 
because we do not restrict damped Lya absorption to ha- 
los, we do not assume any single relation between DLA 
cross section and halo mass for analysis. There is a rela- 
tion between DLA cross-section and associated line width 
which in some cases is indeed caused by the velocity dis- 
persion within individual galaxies. Our findings are as 
follows. 

• Resolution: As we increase grid resolution, our 
column densities demonstrate a trend in conver- 
gence towards the observed distributions. How- 
ever, we cannot reproduce the high-end tail of the 
velocity width distribution. Although it is possi- 
ble that the resolution effects and limitation in box 
sizes are responsible for this shortcoming, feedback 
from star formation which we have not included 
into the current models is likely to produce expand- 
ing shells and therefore higher velocity widths on 
average. The column density distributions alone in- 
dicate that the minimum required grid resolution 
has to be better than ^ 500 pc in comoving coordi- 
nates, or 100 pc in physical space. The minimum 
mass resolution is ^ 10^ corresponding to at 
least a hundred DM particles per halo which is still 
capable of being self-shielded from the UVB. 

• Filamentary DLAs: The nature of DLA ab- 
sorption is a function of halo environment, at 
any given redshift. Although all self-shielded ha- 
los can give rise to DLAs, not all DLAs can be 
associated with individual halos. We find that 
in massive cluster environments at z = 3 tidal 
tails and quasi-filamentary structures often appear 
to be self-shielded against the UVB. The typi- 
cal column densities of filamentary DLAs are be- 
low 10^^ cm^^, and physical widths are of order 
5 — 10 kpc. Since we have not achieved numeri- 
cal convergence, we have not attempted a compre- 
hensive, redshift-dependent statistical study to dis- 
tinguish halo DLAs from the ones caused by tidal 
tails, and both of these from DLAs caused by pri- 
mordial gas still accreting onto galaxies. 

Provided that filamentary structures are still self- 
shielded at higher grid resolution, it seems plausible 
that going to higher resolution would only increase 
the filamentary contribution to DLA absorption, 
as the 3D density peaks (halos) would contribute 
progressively less to the total HI column density 
than the 2D density peaks (filaments) as baryons 
cool and collapse on finer grids. 



• Photoionization of low-mass galaxies: The 

least massive single galaxy associated with DLAs 
in our models has the mass 3.6 x 10^ M0. However, 
we see a significant population of halos with HI in 
absorption only above the mass ~ 7 x 10*" Mq below 
which just a few dwarf galaxies retained any suffi- 
cient amount of neutral gas by z = 3. This thresh- 
old is slightly lower than the pre viously reported 
value s of 10® Mq to few 10® Mq l|Nagamine et alJ 
120041) contributing to larger DLA counts in less 
massive halo environments. Part of this discrep- 
ancy can be explained by spatial variations in the 
UVB inside the galaxy cluster, and it is possible 
that supernova feedback which is not included into 
our models could also have a role removing neutral 
gas from low-mass galaxies. 

• Self-shielding model: Calculations with self- 
shielding above density 10^^ cm~'^ produce roughly 
the same total number of DLAs per unit absorp- 
tion distance as the run with full radiative transfer. 
However, the absorption line width distributions in 
the two models are markedly different due to vari- 
ations in ionizational profiles. 

We intend this work to be the first in a series of pa- 
pers exploring DLA properties with our new cosmological 
simulations with high-angular resolution radiative trans- 
fer on adaptively refined meshes. A number of topics are 
open to further exploration. 

1) At the moment it is not clear to what extent 
DLAs produced with full coupled RHD simulations of 
galaxy formation would be different from the ones we 
get with our radiative-chemical postprocessing approach. 
Inclusion of radiative terms into hydrodynamical equa- 
tions might be important for both accessing the ef- 
fect of the external U V irradiation of primordial galax- 
ies ( I liev et ail l2005tl and studying fee dback from SF 
I^Wha,TeTiTr aLll2004HO'Shea et a].ll2nf)5j) . 

2) We have not included SF into our current mod- 
els. Its mechanical and radiative feedback could change 
the ionization structure and kinematics of gas in high- 
redshift galaxies. For example, winds from SF regions 
and/or SNe could create denser shells around low-density 
bubbles. If dense enough, these shells might contribute 
to an overall increase in DLA effective cross-sections 
(Schavc 2001). Of course, there is a competing effect 
if feedback destroys ne utral gas in halos as observed in 
iNagamine et al.l l|2004D . 

3) Besides an obvious merit of providing more statis- 
tics, going to larger simulation volumes of few tens Mpc 
would allow us to test the effect of large-scale fluctuations 
in the UVB on absorption properties of high-redshift 
galaxies. The average separation between quasars at 
z = 3 — 4 is more than an order of magnitude greater 
than our current largest simulation volume. Similar to 
the proximity effect observed in quasar spectra, very few 
halos would produce a DLA in a cluster that hosts a 
quasar. As we move further away from a quasar, the 
mean background drops. It is also attenuated in large 
clusters of galaxies due to the overall higher IGM column 
densities. If the typical UVB in the vicinity of DLAs is 
50-75% lower than the cosmic average, it could yield dra- 
matically different DLA counts. In addition, larger sim- 
ulation volumes would produce more massive clusters of 
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galaxies, with more violent interactions and more promi- 
nent tidal tails which would also have a direct impact on 
DLA observables. 

4) To restrict our parameter space, we omitted an in- 
teresting topic of fine-tuning the initial power spectrum 
(particularly on galactic scales) to match the observed 
line density of DLAs which we are planning to do in the 
future. 

In order to derive a single realistic population of hydro- 
gen absorbers from Lya forest clouds to DLAs, ideally 
one would like to build a sub-kpc resolution model with 
both stellar feedback and UVB radiative transfer in a 
large (50-80 Mpc) volume. To achieve the required mass 
and grid resolution, one would need to use grids of sizes 
2048L6 with 2048^ DM particles. An N-body model of 
this size (but in a much larger volume) has been recently 
computed by the Virgo Consortium ~ the "Millennium 



Simulation" ijSpringel et al.ll2005(l - and can be poten- 
tially postprocessed with our fully threaded transport 
engine following the routine utilized in this paper. 
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